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Abstract 

In  this  paper  we  develop  local  discontinuous  Galerkin  (LDG)  methods  for 
the  fourth-order  nonlinear  Cahn-Hilliard  equation  and  system.  The  energy 
stability  of  the  LDG  methods  is  proved  for  the  general  nonlinear  case.  Nu¬ 
merical  examples  for  the  Cahn-Hilliard  equation  and  the  Cahn-Hilliard  system 
in  one  and  two  dimensions  are  presented  and  the  numerical  results  illustrate 
the  accuracy  and  capability  of  the  methods. 
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1  Introduction 


In  this  paper,  we  consider  numerical  methods  in  a  bounded  domain  hi  G  (d  <  3) 
for  the  Cahn-Hilliard  equation 

ut  =  V-  (&(u)V(-7Au  +  ^'(u))),  (1.1) 

and  the  Cahn-Hilliard  system 

\ut  =  V  •  (B(«)Vw), 

<  (1-2) 

I  cj  =  — yAw  +  Zf'F('u), 

where  (-ix) and  7  is  a  positive  constant.  Here  b{u)  is  the  non- negative 
diffusion  mobility  and  T  (u)  is  the  homogeneous  free  energy  density  for  the  scalar 
case  (1.1).  For  the  system  case  (1.2),  B{u )  is  the  symmetric  positive  semi-definite 
mobility  matrix  and  \h(w)  is  the  homogeneous  free  energy  density. 

We  develop  a  class  of  local  discontinuous  Galerkin  (LDG)  methods  for  these 
nonlinear  equations.  Our  proposed  scheme  is  high  order  accurate,  nonlinear  stable 
and  flexible  for  arbitrary  h  and  p  adaptivity.  The  proof  of  the  energy  stability  of 
the  scheme  is  given  for  the  general  nonlinear  solutions. 

The  Cahn-Hilliard  equation  was  originally  propose  by  Cahn  and  Hilliard  [8]  to 
study  the  phase  separation  in  binary  alloys.  The  Cahn-Hilliard  system  was  pro¬ 
posed  by  Morral  and  Cahn  [27]  to  model  three-component  alloys.  When  a  single 
homogeneous  system  composed  of  two  or  three  components  at  high  temperature  is 
rapidly  cooled  to  a  temperature  9  below  the  critical  temperature  9C,  the  phase  sep¬ 
aration  happens.  The  Cahn-Hilliard  equations  have  been  adopted  to  model  many 
other  physical  situations,  e.g.  interface  dynamics  in  multi-phase  fluids. 

There  have  been  many  algorithms  developed  and  simulations  performed  for  the 
Cahn-Hilliard  equations,  using  finite  element  methods  [2,  3,  4,  6,  7,  15,  16,  17,  20], 
discontinuous  Galerkin  methods  [9,  21,  31],  multi-grid  method  [23,  24,  25]  and  finite 
difference  methods  [19,  22,  30]. 

Here  we  should  mention  the  difference  between  our  LDG  method  and  the  dis¬ 
continuous  Galerkin  methods  in  [9,  21,  31].  The  discontinuous  Galerkin  method 
considered  in  [9]  refers  to  a  discontinuous  Galerkin  discretization  in  time,  hence  is 
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different  from  our  approach  of  using  a  local  discontinuous  Galerkin  discretization  for 
the  spatial  variables.  The  discontinuous  Galerkin  method  in  [31]  used  the  standard 
C°  finite  element  shape  functions  instead  of  the  discontinuous  basis  functions  in  our 
LDG  method  which  are  allowed  to  be  completely  discontinuous  across  element  in¬ 
terfaces.  In  [21],  a  discontinuous  Galerkin  method  which  is  in  the  DG  family  known 
as  the  interior  penalty  method  [1]  was  developed  for  the  constant  mobility  case  (i.e. 
b[u)  =  constant ).  Stability  was  proved  in  [21,  31],  but  only  for  the  constant  mobility 
case.  Our  LDG  method  does  not  contain  mesh  dependent  stabilization  coefficients 
as  in  [21],  Moreover,  we  prove  stability  for  quite  general  nonlinear  cases,  for  any 
orders  of  accuracy  on  arbitrary  triangulations  in  any  space  dimension. 

The  discontinuous  Galerkin  (DG)  method  is  a  class  of  finite  element  methods, 
using  discontinuous,  piecewise  polynomials  as  the  solution  and  the  test  space.  It  was 
first  designed  as  a  method  for  solving  hyperbolic  conservation  laws  containing  only 
first  order  spatial  derivatives,  e.g.  Reed  and  Hill  [28]  for  solving  linear  equations, 
and  Cockburn  et  al.  [12,  11,  10,  13]  for  solving  nonlinear  equations.  It  is  difficult  to 
apply  the  DG  method  directly  to  the  equations  with  higher  order  derivatives.  The 
idea  of  the  LDG  method  is  to  rewrite  the  equations  with  higher  order  derivatives  into 
a  first  order  system,  then  apply  the  discontinuous  Galerkin  method  on  the  system. 
The  design  of  the  numerical  fluxes  is  the  key  ingredient  to  ensure  stability. 

The  first  LDG  method  was  constructed  by  Cockburn  and  Shu  in  [14]  for  solving 
nonlinear  convection  diffusion  equations  containing  second  order  spatial  derivatives. 
Their  work  was  motivated  by  the  successful  numerical  experiments  of  Bassi  and 
Rebay  [5]  for  the  compressible  Navier-Stokes  equations.  Yan  and  Shu  developed 
an  LDG  method  for  a  general  KdV  type  equation  (containing  third  order  spatial 
derivatives)  in  [36],  and  they  generalized  the  LDG  method  to  PDEs  with  fourth  and 
fifth  order  spatial  derivatives  in  [37].  Levy,  Shu  and  Yan  [26]  developed  LDG  meth¬ 
ods  for  nonlinear  dispersive  equations  that  have  compactly  supported  traveling  wave 
solutions,  the  so-called  “compactons” .  More  recently,  Xu  and  Shu  [32,  33,  34,  35] 
further  developed  the  LDG  method  to  solve  many  nonlinear  wave  equations  with 
higher  order  derivatives,  including  the  general  KdV-Burgers  type  equations,  the 
general  fifth-order  KdV  type  equations,  the  fully  nonlinear  K (n,  n,  n )  equations, 
the  generalized  nonlinear  Schrodinger  equations,  the  coupled  nonlinear  Schrodinger 
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equations,  the  Kuramoto-Sivashinsky  equations,  the  Ito-type  coupled  KdV  equa¬ 
tions,  the  Kadomtsev-Petviashvili  equation,  and  the  Zakharov-Kuznetsov  equation. 
A  common  feature  of  these  LDG  methods  is  that  stability  can  be  proved  for  quite 
general  nonlinear  cases.  DG  and  LDG  methods  also  have  several  attractive  prop¬ 
erties,  such  as  their  flexibility  for  arbitrary  h  and  p  adaptivity  and  their  excellent 
parallel  efficiency. 

The  paper  is  organized  as  follows.  In  Section  2,  we  present  and  analyze  the  local 
discontinuous  Galerkin  methods  for  the  Cahn-Hilliard  system.  In  Section  2.1,  we 
review  the  properties  of  the  Cahn-Hilliard  equation  and  the  Cahn-Hilliard  system. 
In  Section  2.2,  we  present  the  local  discontinuous  Galerkin  methods  for  the  Cahn- 
Hilliard  system.  We  prove  a  theoretical  result  of  the  energy  stability  for  the  nonlinear 
case.  Section  3  contains  numerical  results  for  the  nonlinear  problems  which  include 
the  Cahn-Hilliard  equation  and  the  Cahn-Hilliard  system  for  one-dimensional  and 
two-dimensional  cases.  The  numerical  results  demonstrate  the  accuracy  and  capa¬ 
bility  of  the  methods.  Concluding  remarks  are  given  in  Section  4. 

2  The  LDG  method  for  the  Cahn-Hilliard  system 

2.1  Properties  of  the  Cahn-Hilliard  system 

We  consider  the  model  for  phase  separation  of  a  multi-component  alloy  with  N  >  2 
components  in  bounded  domain  12  G  Wl  (d  <  3).  The  system  of  nonlinear  diffusion 
equations  is  given  by 


ut  =  V  •  (B(m)Voj),  (2.1a) 

uj  = —'jA  u  +  D^(u),  (2.1b) 

^  =  =  011(9 ^  (2-lc) 

u(x,  0)  =  u0(x).  (2. Id) 


Here  x  =  (aq,---  ,Xd),  u,u>  G  ( L2(Q))N ,  {DT(w)};  =  9 ,  d 12  is  the  boundary 
of  12  and  v  is  the  normal  vector  to  <912.  B(u )  is  the  N  x  N  symmetric  positive 
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semi-definite  mobility  matrix  and  has  the  form 


N 


—  Bnp(tl)  I —  ^n(^n)  j  &np  J  ^  ^  bq  i^q)  J  fcp('Up) 

,g=i 


(2.2) 


where  5np  is  the  Kronecker  delta. 

For  n  =  (Vu  -  ■  ■  ,VN),€  =  (6,  ■  ■  ■  ,  &v)  e  ( L2(fl))N  and  S  =  (s^  •  •  •  ,  sN)T,  P  = 
(Pi,--  -  ,2hv)T  with  e  (L2(Q))d,  l  =  1,  •  •  •  ,  TV,  we  set 

M;  =  Vi,  1 7777}  {V?7}z  =  V-/7;,  {Ar7}i  =  A^,  77  ■  £  = 

L  J  ;  2=1 

N 

V  •  S  =  (v  ■  sir  ■  ■  ,V  ■  SN)T,  V  •  5  =  (V  •  Si,  •  •  •  ,  V  •  SAr)T,  S.P  =  y^si.p, 

2=1 

The  concentration  of  the  Ith  component  of  the  alloy  is  denoted  by  Ui  and  so  the 
constraints 

N 

(a)  0<u,  <1,  ( b )  XM  =  1  (2-3) 

2=1 

are  satisfied. 

The  chemical  potential  u:  can  be  dehned  as  the  variational  derivative  of  the 
Ginzburg- Landau  free  energy 


£(u)  ■=  /(-IVu|2  +  ■a(u))dx, 


(2.4) 


i.e.  ui  =  g.  The  gradient  energy  coefficient  7  >  0  and 


'F(u)  :=  vh1(w)  —  ^u1  Au 


(2.5) 


is  the  homogeneous  free  energy  density.  Here,  A  is  a  constant  N  x  N  symmetric 
matrix  taking  into  account  the  interaction  between  different  components.  The  term 
T1(7i)  represents  the  entropy  of  the  system  and  is  usually  taken  to  be  of  the  form 

N 

:=  dyujlnui,  (2.6) 

2=1 

with  the  absolute  temperature  6  >  0.  In  the  deep  quench  limit  9  — »  0,  we  take 


*i(u) 


{0  when  u  satisfies  the  constraints  (2.3), 
00  otherwise. 


(2.7) 
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From  the  boundary  conditions  (2.1c)  we  have 

-f-  f  udx  =  0,  -f- £(u )  <  0.  (2.8) 

dt  Jn  ’  dt  v  v  ' 

Hence,  the  total  mass  of  each  component  is  conserved  and  the  free  energy  £  decays 
for  the  system. 

Remark  2.1.  The  scalar  Cahn-Hilliard  equation  (1.1)  is  a  special  case  of  the  Cahn- 
Hilliard  system  (2.1). 

In  the  case  N  =  2,  assuming  that  An  =  A2 2,  Bn  =  B22,  dehning  u  :=  u2  —  U\, 
to  u2  —  u>i,  b(u)  =  B22  —  B\2  and  6C  =  A22  —  A12,  we  obtain  that  (u,  u)  satishes 
the  equation 

ut  —  V  •  (b{u)Vu)  =0,  c o  —  — 7A u  +  ^(w),  (2.9) 

i.e. 

ut  —  V  ■  ^6(rt)V(— yAu  +  \H'(m))^  =  0,  (2.10) 

with  the  homogeneous  free  energy 

*(«)  =  |  ((1  +  u)  ln(i±b  +  (1  -  u)  ln(  1-^)1  +  |(1  -  u2).  (2.11) 

This  is  the  Cahn-Hilliard  equation  with  a  logarithmic  free  energy  which  satishes  the 
constraint  \u\  <  1. 

We  can  also  define  u  :=  u2  and  w  :=  ,  then  we  obtain  the  same  equation 

(2.10)  with  another  homogeneous  free  energy 

9  9 

T(u)  =  -  (u  Inu  +  (1  —  u)  ln(l  —  u))  +  ~^u(l  —  u),  (2-12) 

which  satishes  the  constraint  0  <  u  <  1. 

The  Ginzburg-Landau  free  energy  of  the  equation  (2.10) 

£{u):=  I  ^|Vu|2  +  T(u))d;r  (2.13) 

also  satishes 

A(„.)  <  0.  (2.14) 
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2.2  The  LDG  method  for  the  Cahn-Hilliard  system 

In  this  section,  we  consider  the  local  discontinuous  Galerkin  method  for  the  Cahn- 
Hilliard  system  (2.1)  with  N  components  in  12  G  Wl  with  d  <  3.  Although  we  do  not 
address  the  numerical  results  in  three  dimensions  in  this  paper,  the  LDG  methods 
and  the  energy  stability  results  of  this  paper  are  valid  for  all  d  <  3. 

2.2.1  Notation 

Let  Th  denote  a  tessellation  of  12  with  shape-regular  elements  K.  Let  T  denote  the 
union  of  the  boundary  faces  of  elements  K  e  %  ,  i.e.  T  =  UKerhdK,  and  r0  =  T\<9f2. 

In  order  to  describe  the  flux  functions  we  need  to  introduce  some  notations.  Let 
e  be  a  face  shared  by  the  “left”  and  “right”  elements  Kr  and  Kr.  For  our  purpose 
“left”  and  “right”  can  be  uniquely  defined  for  each  face  according  to  any  fixed  rule, 
see,  e.g.  [36]  for  more  details  of  such  a  definition.  Define  the  normal  vectors  Ur 
and  Ur  on  e  pointing  exterior  to  Kr  and  Kr,  respectively.  If  -0  is  a  function  on  Kr 
and  Kr,  but  possibly  discontinuous  across  e,  let  'i/jr  denote  {'ip\KL)\e  and  'ipR  denote 
{Mi<R)\ei  the  left  and  right  trace,  respectively. 

Let  VP(K)  be  the  space  of  polynomials  of  degree  at  most  p  >  0  on  K  e  Th.  The 
finite  element  spaces  are  denoted  by 

Vf={V:  <fi\Ke(Vr(K))N,  VKe%}, 

>=«>1,---,4v)T:  [E(W)'  l  —  1-  ■  ■  N,  VK 

Note  that  functions  in  V)'Y  and  E0  are  allowed  to  be  completely  discontinuous  across 
element  interfaces. 

2.2.2  The  LDG  methods 

To  define  the  local  discontinuous  Galerkin  method,  we  rewrite  (2.1)  as  a  first  order 
system: 


e 

<1 

• 

(2.15a) 

S  =  B(u)P, 

(2.15b) 

P  =  V(— q  +  r), 

(2.15c) 
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q  =  7V • W, 
XV  =  Vw, 
r  =  DV(u), 


(2.15d) 

(2.15e) 

(2.15f) 


where  we  use  the  notations  which  are  defined  in  Section  2.1. 

To  simplify  the  notation,  we  still  use  u ,  S,  P,  q,  XV  and  r  to  denote  the 
numerical  solution.  The  local  discontinuous  Galerkin  method  to  solve  the  system 
(2.15)  is  as  follows:  Find  u,  q,  r  e  V);v  and  S,  P,  W  E  £^,  such  that,  for  all  test 
functions  p,  <p,  $,  G  and  ©,  <E>,  T  e  £^, 


/  ut  ■  pdK  —  —  I  S  •  V  pdK  +  v  •  S  ■  pds, 

I K  J  K  JdK 

(  S  •  QdK  =  [  (B(u)P)  •  QdK, 


IK 


IK 


P  •  <S>dK  = 


IK 


f  (r  —  q)  ■  (V  •  <X>)dK  +  f  (r  —  q)  ■  (u  •  $>)ds, 

J  K  JdK 


'K 


q  ■ « pdK  =  —7  XV  •  VipdK  +  7  /  v  •  XV  ■  <pds, 
J  K  JdK 


IK 


W  •  TdK  =  -  u  ■  (V  •  T )dK  +  u  (v  T )ds, 

J  K  JdK 


[  r  ■  £dK  —  [  (DV(u))-£dK. 


(2.16a) 
(2.16b) 
(2.16c) 
(2.16d) 
(2.16e) 
(2. 16f ) 


JK  JK 

The  “hat”  terms  in  (2.16a)-(2.16f)  in  the  cell  boundary  terms  from  integration  by 
parts  are  the  so-called  “numerical  fluxes” ,  which  are  functions  defined  on  the  edges 
and  should  be  designed  based  on  different  guiding  principles  for  different  PDEs  to 
ensure  stability. 

It  turns  out  that  we  can  take  the  simple  choices  such  that 


51,  =  5, 


Q\e  =  qR, 


=  rR, 


xv\e  =  XV 1 


u  e  =  UR. 


(2.17) 


We  remark  that  the  choice  for  the  fluxes  (2.17)  is  not  unique.  In  fact  the  crucial 
part  is  taking  S  and  q,  r  from  opposite  sides  and  XV  and  u  from  opposite  sides. 

Remark  2.2.  For  the  scalar  Cahn-Hilliard  equation 

ut  =  V  ■  [b{u)V[— 7 Aw  +  'F'('u))^, 
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(2.18) 


the  LDG  scheme  becomes:  Find  u,  q,  r  G  V f  and  s,  p,  w  G  such  that,  for  all 
test  functions  p,  99,  £  G  V),1  and  Tj,  (j),  -0  G  E^ 


/  utpdK  =  —  s  ■  V pdK  +  /  uT~spds, 

J  K  J  K  JdK 

/  s  ■  rjdK  =  /  6(u)p  ■  rjd.K , 

J  I<  J  K 

/  p  ■  —  —  I  (r  —  g)V  ■  4>dK  +  (r  — 

•7/c  Ja'  iax 

/  qpd.K  =  —7  /  in  •  VpdK  +  7  /  C^wpds, 

J  K  J K  JdK 

/  in  ■  t^dK  —  —  I  mV  ■  ifdK  +  /  mi/  •  tfds, 

/  A'  J  K  JdK 

[  r£dK  —  I  ^'{u)£dK. 

J  K  J  K 


(2.19a) 

(2.19b) 

(2.19c) 

(2.19d) 

(2.19e) 

(2.19f) 


The  numerical  fluxes  are 


«|e  =  «L,  51  e  =  qR,  r\e  =  r  Ri  w\e  =  WL,  u\e  =  Ur. 


(2.20) 


2.2.3  Energy  stability 

We  will  prove  the  theoretical  results  of  the  energy  stability  for  the  general  nonlinear 
system  case  with  the  choice  of  the  fluxes  in  the  previous  section. 


Proposition  2.1.  (Energy  stability)  The  solution  to  the  schemes  (2.16)  and  (2.17) 
with  the  boundary  conditions  (2.1c)  satisfies  the  energy  stability 

T  [  (eW*  W  +  V(u))dx  <  0. 
dt  Jn\ 2  / 

Proof.  Choosing  the  test  function  ^  =  —  ut  G  Vff  in  (2. 16f) ,  we  obtain 

-  j  r  ■  utdK  =  —  [  (Dm(u))-utdK.  (2.21) 

J  K  Jk 

For  the  equation  (2.16e),  we  first  take  the  time  derivative,  then  choose  the  test 
function  T  =  7 W  G  E^.  We  obtain 

7  f  Wt  •  WdK  =  -7  /  ut  ■  (V  •  W)dK  +  7  f  ut  ■  (v  •  W)ds.  (2.22) 
•7  A'  JA'  JdK 

For  (2.16a),  (2.16b),  (2.16c)  and  (2.16d),  we  take  the  test  functions 
p  =  r  -  q,  ©  =  -P,  $  =  5,  <p  =  ut. 
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Then  we  have 


/  ut  ■ (r  —  q)dK  =  —  (  S  •  (V(r  —  q))dK  +  f  is  •  S  •  (r  —  qr)ds, 

./AT  JK  JdK 

-  I  S  •  PdK  =  -  I  (B(u)P)  •  PdP, 

Ja'  Ja' 

/  P  •  SdP  =  -  [  (r  ~  g)  •  (V  •  5)dP  +  [  (r-q)-{v  S)ds, . 

J  K  J K  JdK 

/  q  ■  utdK  =  —7  W  •  (’ Vut)dK  +  7  /  is  •  W  ■  utds. 

J  K  Jk  JdK 

Summing  up  the  equations  (2.21)-(2.26),  we  obtain 

J  (jW  •Wt  +  (DV(u))  ■ut'j  +  J  (P(w)P)  •  PdK 

=  -7  [  ut(V  mW)dK  -7  [  W»{Vut)dK 
Jk  Jk 

+  7  /  ut  ■  (is  •  W)ds  +  7  /  is  •  W  ■  utds 

JdK  JdK 

-  [  S  •  (V(r  -  q))dK  -  [  (r  -  q)  ■  (V  •  S)dK 


(2.23) 

(2.24) 

(2.25) 

(2.26) 


'  A' 


'  AT 


+  /  •  5  •  (r  —  g)ds  +  /  (r  —  q)  ■  (is  •  S)ds, 

JdK  JdK 

=  -7  /  (P*  W)-'Ujcis  +  7  /  tij  •  (is  •  W)ds  +  7  /  is  •  W  ■  utds 

JdK  JdK  JdK 

—  f  (r  —  q)  ■  (is  •  S)ds  +  f  is  •  S  ■  (r  —  q)ds  +  f  (r  —  q)  ■  (is  •  S)ds. 

JdK  JdK  JdK 

Summing  up  over  Jl ,  with  the  numerical  fluxes  (2.17)  and  the  boundary  conditions 
(2.1c),  we  get 

J  (jW  •Wt  +  ty(u)ty,x  +  J  (B(u)P)  •  Pdx  =  0. 

Because  B(u)  is  semi-positive,  we  have  the  energy  stability 

T  f  (t:W  +  ^(u))dx  <  0. 
dt  Jq\2  J 


□ 


Remark  2.3.  Proposition  2.1  is  also  true  for  the  LDG  scheme  (2.19)  and  (2.20)  for 
the  scalar  Cahn-Hilliard  equation  (2.18).  The  proof  goes  along  the  same  line  and  is 
simpler.  We  thus  omit  the  details. 
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3  Numerical  results 


In  this  section  we  perform  numerical  experiments  of  the  local  discontinuous  Galerkin 
method  applied  to  the  Cahn-Hilliard  equation  and  system.  Time  discretization  is  by 
the  third  order  TVD  Runge-Kutta  method  [29].  We  have  chosen  At  suitably  small 
so  that  spatial  errors  dominate  in  the  numerical  results.  This  is  not  the  most  efficient 
method  for  the  time  discretization  to  our  LDG  scheme.  However,  we  will  not  address 
the  issue  of  time  discretization  efficiency  in  this  paper.  All  the  computations  were 
performed  in  double  precision.  We  have  verified  with  the  aid  of  successive  mesh 
refinements,  that  in  all  cases,  the  results  shown  are  numerically  convergent. 

3.1  Numerical  results  for  the  Cahn-Hilliard  equation 

3.1.1  One  space  dimension 

In  this  section,  we  give  the  numerical  test  results  for  the  one-dimensional  Cahn- 
Hilliard  equation. 


Example  3.1. 

We  consider 


ut  =  ~(b(u)(^uxxx  -  (T'O))^  (3.1) 

with  T(w)  =  |(1  —  u2),  b{u )  =  1  —  u2  and  7  =  0.01  in  H  =  (0,1).  The  initial 
condition  is 


u0{x) 


cos 


vT 


-  1, 


if  \x  _  il  <  AH 
11  r  2I  —  2  ’ 

otherwise. 


-1, 


(3.2) 


The  boundary  conditions  are  taken  as 


v>x  b{v?)uxxx  0  (3.3) 

at  both  ends.  We  note  that  uq{x)  is  in  W1( fl)  and  not  in  H2(Q).  Elliot  and  Garcke 
[18]  proved  existence  of  a  solution  with  the  property  that  u  G  L2(0,  T;  H2(Q))  for 
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Table  3.1:  Accuracy  test  for  the  Cahn-Hilliard  equation  (3.1)  with  the  stationary 
solution  (3.4).  Uniform  meshes  with  J  cells  at  time  t  =  0.1. 


J 

L°°  error 

order 

L 2  error 

order 

10 

1.85E-01 

- 

6.94E-02 

- 

pO 

20 

1.44E-01 

0.37 

4.44E-02 

0.64 

40 

6.83E-02 

1.07 

2.08E-02 

1.09 

80 

2.97E-02 

1.19 

8.67E-03 

1.26 

10 

7.55E-02 

- 

2.42E-02 

- 

p 1 

20 

1.45E-02 

2.38 

3.86E-03 

2.64 

40 

4.06E-03 

1.83 

8.11E-04 

2.25 

80 

9.07E-04 

2.16 

1.96E-04 

2.04 

arbitrary  initial  data  Uq  G  Hl(Vt).  Onr  numerical  tests  verify  their  conclusion  that 
the  numerical  solution  appears  to  spread  to  the  stationary  C'1([0, 1])  solution: 


^ steady  (A) 


1  +  cos(^ 


-1, 


-1,  if  \x-  ||  <7T^ 
otherwise. 


(3.4) 


The  L 2  and  L°°  errors  and  the  numerical  orders  of  accuracy  for  the  stationary 
solution  u steady  at  time  t  —  0.1  with  uniform  meshes  in  [0, 1]  are  contained  in  Table 
3.1.  In  Fig.  3.1,  we  show  the  numerical  results  at  t  =  0.1  using  P 1  elements  on 
the  uniform  mesh  with  80  cells.  With  fewer  cells,  onr  scheme  gets  the  same  results 
comparing  the  numerical  calculations  performed  by  Barrett  el  al.  [3]. 


Example  3.2. 

We  consider  the  Cahn-Hilliard  equation  (3.1)  with  b(u)  =  1  or  b{u )  =  1  —  u2  and 
7  =  10~3  in  fl  =  (0, 1).  We  take  the  free  energy 


T(u)  =  °- 


(1  +  u)  ln(— -7— -)  +  (1  -  u)  ln(— — — ) 


+  2^ 


(3.5) 
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Figure  3.1:  The  numerical  solution  of  the  equation  (3.1)  with  the  initial  condition 
(3.2)  and  the  boundary  conditions  (3.3)  at  t  —  0.1  using  P 1  elements  on  the  uniform 
mesh  with  80  cells. 


with  9  =  0  (the  deep  quench  limit)  or  0.3.  The  initial  condition  is 


u0(x)  =  { 


20(|  -x), 

-20k -SI, 

-1, 


if  0  <  x  <  l  —  S, 


if  It  —  —  I  ^  — — 

3  I  —  20’ 

if  411  <  1_ 
11  I**'  50  I  —  20’ 


50  I 

otherwise. 


(3.6) 


The  boundary  conditions  are  (3.3). 

We  use  P 1  element  and  a  uniform  mesh  with  80  cells.  The  results  include  both 
6  =  0  (the  deep  quench  limit)  and  9  =  0.3  for  constant  and  degenerate  mobility 
b(u)  =  1  or  b(u)  =  1  —  u2.  The  simulations  are  stopped  when  the  obtained  profiles  do 
not  change  for  a  long  time.  The  numerical  results  compare  very  well  with  numerical 
calculations  performed  by  Barrett  el  al.  [3].  From  the  numerical  results  in  Fig.  3.2, 
we  have  the  following  observation: 


•  For  the  constant  mobility  b(u)  =  1,  the  “bump”  is  swept  away  quickly.  This 
is  due  to  the  fact  that  mobility  is  positive  in  the  pure  phases. 


•  For  the  degenerate  mobility  b(u)  =  1  —  u2  with  logarithmic  free  energy  (3.5), 
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0=0,  b(u)=1 


0=0,  b(u)=1-u2 


0  0.2  0.4  0.6  0.8  1  0  0.2  0.4  0.6  0.8  1 


Figure  3.2:  The  solution  of  the  equation  (3.1)  with  the  initial  condition  (3.6)  and 
the  boundary  conditions  (3.3)  at  different  time  T  with  P 1  element  on  the  uniform 
mesh  with  80  cells. 

the  time  scale  of  the  diffusion  is  greatly  increased. 

•  For  the  degenerate  mobility  b(u)  =  1  —  u2  and  the  quench  limit  free  energy, 
the  “bump”  does  not  lose  mass.  As  6  goes  to  zero,  the  minima  of  the  free 
energy  T(w)  in  (3.5)  converge  to  u  —  ±1  (see  Fig.  3.3).  This  implies  that  the 
diffusion  through  the  bulk  becomes  smaller  for  lower  temperature. 

3.1.2  Two  space  dimensions 

In  this  section,  we  present  the  numerical  results  for  the  two-dimensional  Cahn- 
Hilliard  equation. 
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Figure  3.3:  The  free  energy  T(u)  in  (3.5)  with  9  =  0.3,  0.5,  0.6. 


Example  3.3. 

We  consider  the  Cahn-Hilliard  equation 

ut  =  V  •  (&(u)V(-7Au  +  tf'(u)))  (3.7) 

with 


T(u)  =  600(u Inn  +  (1  —  u)  ln(l  —  u))  +  1800w(l  —  u),  b{u)  =  1,  7  =  1. 


The  initial  condition  is 


«o(x)  = 


0.71  x  G  Oi , 

0.69  x  G  1^2, 


where  the  square  domain 


(3.8) 


n=  (-0.5, 0.5)  x  (-0.5, 0.5),  S7i  =  (—0.2, 0.2)  x  (—0.2, 0.2),  Sl2  =  Sl-fii. 

The  boundary  conditions  are 

Qu 

—  =  5(m)Vo;  •  v  =  0,  on  d fl.  (3.9) 

We  use  the  P°  and  P 1  element  on  the  uniform  meshes  with  40  x  40  and  80  x  80 
cells  respectively.  The  contours  at  t  —  8  x  10~5  are  shown  in  Fig.  3.4.  We  can  see 
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P°,  40x40 


P\  40x40 


P°,  80x80 


P1,  80x80 


Figure  3.4:  The  contours  of  u(x,t)  for  the  equation  (3.7)  with  the  initial  condition 
(3.8)  and  the  boundary  conditions  (3.9)  when  t  —  8  x  10“5.  P°  and  P 1  elements  on 
the  uniform  mesh  with  40  x  40  and  80  x  80  cells. 

that  the  solution  structure  is  well  resolved  even  for  the  coarser  mesh.  The  numerical 
results  compare  very  well  with  the  numerical  calculations  performed  by  Wells  et  al. 

[31]- 

Example  3.4. 

In  the  square  domain  fl  =  (—0.5,  0.5)  x  (—0.5,  0.5),  we  consider  the  Cahn- Hilliard 
equation  (3.7)  with 

T(n)  =  3000(wlnu  +  (1  —  u)  ln(l  —  u))  +  9000w(l  —  u),  b(u)—u(l  —  u),  7  =  1. 
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Figure  3.5:  The  contours  evolution  of  u(x,t)  for  the  equation  (3.7)  at  different  time 
from  a  randomly  perturbed  initial  condition  with  Pl  element  on  the  uniform  mesh 
with  80  x  80  cells. 

The  initial  condition  Uq  is  a  random  perturbation  of  uniform  state  u  =  0.63  with 
a  fluctuation  no  larger  than  0.05.  The  boundary  conditions  are  taken  as  (3.9). 
This  example  is  used  in  Sec.  5.3  in  [31]  (the  initial  condition  is  identical  in  the 
statistical  sense).  We  use  the  P 1  element  on  a  uniform  mesh  with  80  x  80  cells. 
The  concentration  evolution  can  be  categorized  in  two  stages.  The  first  stage  is 
governed  by  spinodal  decomposition  and  phase  separation  (the  first  four  figures  in 
Fig.  3.5).  The  second  stage  is  governed  by  grain  coarsening  (from  t  —  8  x  10-6 
onwards).  Fig.  3.5  shows  statistically  similar  patterns  in  the  numerical  solution  as 
those  in  Wells  et  al.  [31]. 
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3.2  Numerical  results  for  the  Cahn-Hilliard  system 


3.2.1  One  space  dimension 

In  this  section,  we  present  the  numerical  experiment  results  for  the  one-dimensional 
Cahn-Hilliard  system. 

Example  3.5. 

We  consider  a  ternary  system  in  12  =  (0, 1)  by  Blowey  et  ah  [7] 

Ut  T  t  Cr.r.r.r  0CUXX  .B{Z2\lH(li)}a:x  0,  (3.10) 


with 


and 


B 


(  2  _1  _l\ 

3  3  3 

_1  2  _  1 
3  3  3 

_  1  _  1  2 

V  3  3  3  / 


H/(w)  =  d(ui  In  Mi  +  U2  In  u2  +  u3lnu3)  +  0c(uim2  +  W2W3  +  w3wi). 


The  boundary  conditions  are 


uT  =  Burrx  =  0 


(3.11) 


at  both  ends. 

We  first  perform  a  linear  stability  analysis.  We  seek  a  solution  of  the  form 

OO 

Ui(x,  t )  =  rrii  +  c"(2)  cosnnx,  i  —  1,  •  •  •  ,  3, 

n=l 

where  m  =  (roi,m2,ro3)  is  the  mean  concentration  and  |cf(£)|  <C  1.  Note  that 
mi  +  m2  +  m3  =  1  and  c"(t)  +  c%(t)  +  c3 (t)  =  0.  Linearizing  D^>i{u)  about  rrq  and 
substituting  into  (3.10),  we  obtain  the  ordinary  differential  equations 

dcr‘ 


where 


cn(t)  —  (ci ,  c%),  H  = 


dt 


20  1 


n47r47Cn  +  n27r2Hcn  =  0, 


(3.12) 


+ 


i 


3  l  mi  2(1— mi— m2) 


_20  ( _ 1 _ 


29  (  _L _ 1 

3  \  2m2  2(1— mi— m2) 


1 


20  1 


3  y  2mi  2(1— mi— m2)  J  3  y  m2  2(1— mi— m2) 


+ 


-0r. 
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The  solution  of  (3.12)  is  given  by 


\t)  =  e~n4n4'lt 


—  TL^TT^Ut  T 

x  e  c 


(0). 


For  the  growth  of  one  or  more  of  the  components  u \ ,  w2,  a  necessary  condition  is 
that  the  eigenvalues  of  H  is  smaller  than  — jn2.  When  m2  =  mi,  we  have 

dct  H  -  - - —  +  6mi  ~  ~  ml)g2 

3mf(l  —  2mi_)  c 

We  see  from  Fig.  3.6  that  two  curves  =  mi  and  6/9c  =  3 mi  —  6 m\  define  the 
four  regions  where  H  is  positive,  negative  definite  or  indefinite. 


Figure  3.6:  The  positive,  negative  definite  and  indefinite  regions  of  H,  when  m2  = 

mi. 

We  take  6  =  0.2,  9C  =  1  and  7  =  5.0  x  10“3.  The  initial  conditions  are  random 
perturbations  of  the  uniform  state  m  with  the  fluctuation  no  larger  than  0.01.  We 
use  P 1  element  and  a  uniform  mesh  with  80  cells.  The  simulations  are  stopped  when 
the  obtained  profiles  do  not  change  for  a  long  time. 

We  perform  four  experiments  with  initial  data  inside  the  positive,  negative 
definite  and  indefinite  regions  respectively,  by  taking  rrii  =  1/20,3/20,1/3,19/20 
(points  a,  b,  c  and  d  in  Fig.  3.6  respectively) 

•  mi  =  1/20  in  the  positive  definite  region. 

Fig.  3.7  shows  the  time  evolution  of  the  ternary  system  (3.10).  As  expected, 
the  homogeneous  system  is  stable  and  no  phase  separation  happens. 
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Figure  3.7:  The  evolution  of  system  (3.10)  at  different  time  T  with  m\  =  1/20  when 
9  =  0.2  and  9C  =  1. 

•  m i  =  3/20  in  the  indefinite  region. 

Fig.  3.8  shows  the  time  evolution  of  the  ternary  system  (3.10).  Initially  the 
third  phase  113  dominates.  For  some  time  the  evolution  is  in  the  direction  of 
u\  =  U2  with  two-phase  structure. 

•  m  1  =  1/3  in  the  negative  definite  region. 

Fig.  3.9  shows  the  time  evolution  of  the  ternary  system  (3.10).  We  observe 
three  phases  in  the  early  stages  of  the  spinodal  decomposition. 

•  m\  =  19/20  in  the  indefinite  region. 

Fig.  3.10  shows  the  time  evolution  of  the  ternary  system  (3.10).  The  decom¬ 
position  process  is  like  a  binary  alloy.  After  the  quench,  only  U\  and  112  are 
separated  and  there  is  no  spatial  area  where  M3  is  dominant. 


3.2.2  Two  space  dimensions 

In  this  section,  we  present  the  numerical  simulation  results  for  the  two-dimensional 
Cahn-Hilliard  system. 

Example  3.6. 
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Figure  3.8:  The  evolution  of  system  (3.10)  at  different  time  T  with  m\  =  3/20  when 
6  =  0.2  and  6C  =  1. 


We  consider  a  ternary  system  in  =  (0, 1)  x  (0, 1) 

ut  =  V  •  (B(w)Vw), 

u:  =  — jAu  +  .DTi('u)  —  Au 


\ 


where  ^(w)  is  given  by  (2.6)  and 

^ui(w2  +  -u3)  -MiM2  —  »1»3 

B(u)  —  -uiu2  u2(u1  +  u3)  ~u2u3 

l  -uiu3  - u2u3  u3(u1+u2)l 


A  =  -0r 


(3.13) 


^0  1  Is 
1  0  1 
U  1  0 


Figure  3.9:  The  evolution  of  system  (3.10)  at  different  time  T  with  rn \  =  1/3  when 
9  =  0.2  and  9C  =  1. 


We  take  9  =  1200,  9C  =  3600  and  7  =  1.  The  initial  data  is 
' 

(0,  0, 1)T  if  0  <  Xi  <  j|  and  x2  >  0.65  +  cos(87rxi) 

or  <  X\  <  y|  +  and  x2  >  0.65  —  V3(xi  —  j§) 
or  ||  +  ^  <  xi  and  x2  > 


u0(x1,x2)  = 


(0, 1,  Of  if  0  <  x\  <  y|  and  x2  <  0.35  —  cos(87txi) 

or  y§  <  Xi  <  y|  +  and  x2  <  0.35  +  f3f  1  —  j§) 
or  if  +  5^  -  Xl  an<^  ^2  < 

(0, \)T  if  if  +  <  xx  and  x2  = 

( 1.  0.  0)T  otherwise. 


Figure  3.10:  The  evolution  of  system  (3.10)  at  different  time  T  with  m\  =  19/20 
when  9  =  0.2  and  9C  =  1. 

The  boundary  conditions  are 

=  B{u)^-  =  0,  on  dQ.  (3.15) 

ou  ou 

We  show  the  contours  of  ui(x,t),  u2(x,t )  and  u3(x,t)  at  t  =  8  x  10~5  in  Fig.  3.11 
using  P 1  element  on  a  uniform  mesh  with  80  x  80  cells.  As  expected,  the  symmetry 
of  the  initial  data  is  maintained  during  the  evolution.  We  find  that  the  interface  of 
the  two  components  is  “wetted”  by  the  third  component.  This  is  understood  as  the 
energy  required  to  go  directly  from  the  first  to  the  third  component  is  much  greater 
than  that  required  to  go  via  the  intermediate  second  component.  This  phenomenon 
is  known  as  “wetting” . 


23 


u,(t=0) 


u,(t=8E-5) 


Figure  3.11:  The  contours  of  Ui(x,t),  u2{x1t)  and  Us(x,t)  for  the  equation  (3.13) 
with  the  initial  condition  (3.14)  and  the  boundary  conditions  (3.15)  when  t  =  8  x 
1CF5.  P 1  elements  on  the  uniform  mesh  with  80  x  80  cells. 
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4  Conclusion 


We  have  developed  local  discontinuous  Galerkin  methods  to  solve  the  Cahn-Hilliard 
equation  and  the  Cahn-Hilliard  system.  The  energy  stability  is  proven  for  general 
nonlinear  case.  Numerical  examples  for  one- dimensional  and  two  dimensional  cases 
are  given  to  illustrate  the  accuracy  and  capability  of  the  methods.  Although  not 
addressed  in  this  paper,  the  LDG  methods  are  flexible  for  general  geometry,  unstruc¬ 
tured  meshes  and  h-p  adaptivity,  and  have  excellent  parallel  efficiency.  The  LDG 
method  has  a  good  potential  in  solving  the  Cahn-Hilliard  equations  and  similar 
nonlinear  equations  in  mathematical  physics. 
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